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NUMERICAL METHODS IN THE CONTEXT OF COMPARTMENTAL MODELS 

IN EPIDEMIOLOGY 

Peter Kratz^, Etienne Pardoux^ and Brice Samegni Kepgnou^ 


Abstract. We consider compartmental models in epidemiology. For the study of the divergence of the 
stochastic model from its corresponding deterministic limit (i.e., the solution of an ODE) for long time 
horizon, a large deviations principle suggests a thorough numerical analysis of the two models. The 
aim of this paper is to present three such motivated numerical works. We first compute the solution 
of the ODE model by means of a non-standard finite difference scheme. Next we solve a constraint 
optimization problem via discrete-time dynamic programming: this enables us to compute the leading 
term in the large deviations principle of the time of extinction of a given disease. Finally, we apply the 
T-leaping algorithm to the stochastic model in order to simulate its solution efficiently. We illustrate 
these numerical methods by applying them to two examples. 

Resume. On considere des modeles comportementaux en epidemiologie. Afin d’etudier I’ecart en 
temps long entre le modele stochastique et sa limite loi des grands nombres (qui est la solution d’une 
EDO), on se base sur un principe des grandes deviations, qui nous conduit a mener une etude numerique 
des deux modeles, sur trois aspects differents. Tout d’abord, nous calculons une solution approchee de 
I’EDO a I’aide d’une methode numerique dite “non-standard”. Ensuite une resolvons numeriquement 
un probleme de controle sous contrainte, afin de calculer le terme principal des grandes deviations 
du temps de sortie d’une situation endemique. Enfin nous mettons en oeuvre I’algorithme du “r- 
leaping” pour simuler efficacement la solution du systeme stochastique. Nous illustrons ces simulations 
numeriques en les appliquant a deux examples. 


Introdugtion 

It is well-known that deterministic ODE models of population dynamics (such as the evolution of diseases) 
are not appropriate for small population sizes N (i.e., N < 10^, 10®). Recently, it has been shown (see [^) that 
in some cases the ODE models can diverge from the corresponding stochastic models even for larger population 
sizes {N > 10®). 

We consider the “natural” stochastic extensions for many compartmental ODE models in epidemiology: 
individual based Poisson driven models. The ODE models can be recovered from these models as the law 
of large numbers limit as the population size N approaches infinity. Hence, the deterministic model can be 
considered an appropriate first approximation of the stochastic model for large N. While the dynamics of many 
of the underlying deterministic models are well-understood, the analysis of the corresponding stochastic models 
is often difficult. 

We are interested in the long-term behavior of the epidemic process. For many deterministic models, the 
disease will finally either die out or become endemic, depending on the parameter choice and/or the initial 
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sizes of the compartments, i.e., the solution of the ODE converges to an asymptotically stable equilibrium. 
The theory of Large Deviations (also called the large deviations principle (LDP)) allows us to quantify the 
time taken by the stochastic model to diverge significantly from its deterministic limit. In particular, one can 
compute asymptotically, for large N, the expected time which is needed to transfer the epidemic process from 
the domain of attraction of one equilibrium to the one of another equilibrium. For the simplest models, this 
helps, e.g., to answer the question of when an epidemic becomes extinct although it should be endemic according 
to the deterministic approximation. Particularly interesting cases are epidemiological models where the ODE 
has several stable equilibria (usually a disease-free equilibrium and an endemic equilibrium as, e.g., in |^[^). 

The remainder of this article is structured as follows. In section we formulate stochastic and deterministic 
compartmental models and introduce two specific models which we analyze throughout this article by means of 
numerical methods (section [l.l[ ). We furthermore describe how the time of exit from the domain of attraction of 
a stable equilibrium can be computed (section 1.2). This motivates three different numerical projects which we 
outline in section [L3l We address these questions in sections[^-|^as follows. First we apply a non-standard finite 
difference scheme due to in order to compute the solution of the deterministic model numerically. Second, 
in section]^ by solving a control problem via dynamic programming, we compute the leading coefficient in the 
large deviations analysis of the time of exit by the solution of the SDE from the domain of attraction of a 
locally stable equilibrium of the ODE. Finally, we simulate the stochastic process using the so-called r-leaping 
algorithm. 


1. Compartmental models in epidemiology 

1.1. Stochastic and deterministic models 

We consider diseases in large populations of (initially) N individuals which are split up into different groups 
(or compartments) according to the disease status of the individuals. Such groups are, e.g., the group of 
individuals susceptible to the disease, the group of infectious individuals and the group of immune individuals. 
In order to better compare the population dynamics for different values of N, we normalize the compartment 
sizes by dividing by N; instead of considering the number of individuals in the different groups, we hence 
rather consider their “proportions” with respect to the total size of the population (note that for models with 
non-constant population size, these are only initially truly proportions). 

For a closed set A C IR'^, we define the d-dimensional jump process via a set of k jump directions hj G llA 
and respective Poisson rates /3j{z), j = 1, ..., fc; {t) denotes the “proportion” of individuals in compartment 
* = 1,..., d at time t > 0: 

^ "t 

Z^{t)-.= Z^'^{t):=x+^^h,P,n NI3,{Z^{s))ds) (1) 

j=i 

= x+ hMZ^{s))ds + 1 ^ h,M, ( t N/3,{Z^{s))ds) ; 

Jo Jo 

here x G A, {Pj)j=i,...,k are i.i.d. standard Poisson processes and Mj{t) = Pj{t) — t, j = are the 

associated martingales. The jump directions hj and their respective rates are chosen in such a way that 
Z^(t) G A almost surely for all t. 

The dynamics of the deterministic model corresponding to the stochastic process in Q are given by the 
solution of the following ODE: 


Z{t) ■=Z^{t) = x + 



k 

Y,hjl3jiZ{s))ds. 

i=i 


( 2 ) 
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The two models are connected via a Law of Large Numbers (LLN) (see [^): 
continuous, we have 


if the rates Pj are Lipschitz 


lim = Z^{t) a.s. uniformly on compact time intervals [0,T]. (3) 

N—^OO 

Note also that in the case of Lipschitz continuous rates, equation ([^ admits a unique solution on [0,T]. 

Example 1.1. SIS model. We first consider a simple SIS model without demography (S'(t) being the number 
of susceptible individuals and I{t) = N — S{t) the number of infectious individuals at time t) . For /3 > 0, 
we assume that the rate of infections is /3S'(t)/(t)/7Vj^ In a similar way, we assume for 7 > 0 that the rate 
of recoveries is 'yl{t). We ignore immunity due to earlier infection and assume that an infectious individual 
becomes susceptible after recovery. As population size is constant, we can reduce the dimension of the model 
by solely considering the proportion of infectious at time t, that is d = 1 and Zt = It- Using the notation of 
equations 0 and ([^, we have 

A =[0,1], hi = l, /i2 = -l, /3i(z) =/32(1 - z), /32 (z)=7z. 

It is easy to see that the ODE (|^ has a disease free equilibrium x = 0. This equilibrium is asymptotically stable 
if Ro = /3/7 < 10 If i?o > 1, is unstable and there exists a second, endemic equilibrium x* = 1 — 7 //? which 
is asymptotically stable. While in the deterministic model the proportion of infectious individuals converges to 
the endemic equilibrium a;*, the disease goes ultimately extinct in the stochastic model. 

Example 1.2. A model with vaccination. We consider a (deterministic) model with vaccination and 
demography described in and its stochastic counterpart. We illustrate the different transitions in Figure 1; 
here, S and I denote the number of susceptible respectively infectious individuals as before and V denotes the 
number of vaccinated individuals. We assume that individuals are vaccinated at a certain rate but can lose their 
protection again; the vaccine is not perfect but decreases the rate of infection by a factor cr S [ 0 , 1 ]: if cr = 1 , 
the vaccine is useless, if cr = 0 it protects the individuals perfectly. Furthermore, we assume that individuals 
are born susceptible and die (at the same rate) independently of their disease status. This ensures that the 
population size remains constant in the deterministic model. For simplicity, we ensure constant population size 
by synchronizing births and deaths in the stochastic model. Hence, we can again reduce the dimension of the 
models; the first coordinate of the process denotes the proportion of infective individuals, the second coordinate 
denotes the proportion of vaccinated individuals. Using the notations of equations 0 and ([^, we obtain 

A. = {z G 1 R ^|0 ^ zi + Z 2 ^ 1 }, 

hi = (1,0)^, /3i(z) =/?zi(I - zi - Z 2 ), h2 = (I,-l)^, /32iz) = crl3ziZ2, 

ha = (-1,0)^, /33 (z)=7Zi, h4 = (0, I)^, I3i{z) = 9 z2, h5 = (0,-I)^, P^iz) = r){l - zi - Z 2 ), 

he = (-1,0)^, Pti{z) = hr = (0,-1)^, /37(z) = ^Z2- 

The basic reproduction number of this model is given by 

P H + e + ar] ~ P 
” /r + 7 /r + h + ? 7 ’ ° M + t’ 

here, i?o denotes the basic reproduction number without vaccination, i.e., U(0) = 0 and p = Q (see [^). There 
exists a disease-free equilibrium x (xi = 0 ) which is asymptotically stable for Rq < 1 and unstable for Rq > 1 . 

^The reasoning behind this is the following. Each individual (in particular each infected individual) meets other individuals at 
rate a. The probability that the encounter is with a susceptible is S{t)/N. Denote by p the probability that such an event yields 
a new infection. Hence the total rate of new infections is ^S{t)I{t)/N, if /3 = pa. 

^Rq is the so-called basic reproduction number. It denotes the average number of secondary cases infected by one primary case 
during its infectious period, at the start of the epidemic. 
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Figure 1. Schematic representation of the disease transmission for the model with vaccination by [^. 

For Rq < 1 < Rq (and additional assumptions on the parameters, see again), two endemic equilibria x* and x 
{x\,xi > 0) exist: one (say x*) is locally asymptotically stable and one is unstable; the disease-free equilibrium 
is locally asymptotically stable. 


1.2. Large deviations 

In order to quantify the deviation of the stochastic model from the deterministic model, an LDP is desirable. 
For the models we have in mind (cf. Examples [H] and |1.2[ ) , some of the rates Pj : A ^ 1R,+ tend to zero as x 
approaches the boundary of A (this ensures that the process does not leave the domain A). For this reason, 
general results for LDPs for Markov jump processes (see [6||^) do not apply. provides a large deviations 
principle which allows for such diminis hing rates; however, their assumptions do often not apply to the more 
complicated models in epidemiologyj^ provides an LDP on the space D([0, T]; for a large class of 
epidemiological models with compact domain A by generalyzing the results from [^; the rate function It^x of 
the LDP is given as follows}^ 


J/o if is absolutely continuous and (/)(0) = a; 

J- 1^\T J ' \ 1 

cx) else. 


where L denotes the Legendre-Fenchel transform: 


L{x,y) := sup e{p,x,y) 

pGR*^ 


(4) 


for 


e{p,x,y) = {p,y) -'^Pj{x){e^P'^^'> -1). 


i=i 


By d , we have that 

It,x{P) =0 if and only if p solves the ODE Q. (5) 

We interpret It,x{4') as the action required to follow a trajectory </), or in other words the action required to 
deviate from the solution of the ODE. 


^The assumptions in apply to the SIS model (Example 1.1 \ but not to the model with vaccination in Example 


1.2 


^Here, D([Q,T]\ A) denotes the space of cadlag functions on [0, T] with domain A equipped with an appropriate metric such 
that the resul ting space is Polish (see, e.g., [ll| ). 


^See, 


resul ting 

s.g., g 


for an introduction to the theory of large deviations and a definition of the rate function. 
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We are particularly interested in the Freidlin-Wentzell theory (see, e.g., 12 O). First, we want to compute 


the time at which the process leaves the domain of attraction of an asymptotically stable equilibrium 


O, i.e., 


:= inf{t > S ^ \ O}. 

In the case of the SIS model (and Rq > 1) we consider the set O = (0,1]; for the model with vaccination (and 
appropriate parameter choice, in particular i?o < 1 < Ro)j we consider the set 


0 = {z&A\ lim Z^it) = X*} ® 

In other words, we are interested in time of extinction of the diseas^which should eventually become endemic 
according to the deterministic model. 

If A is compact, it can be deduced from the LDP that for all 5 > 0, 


lim pre^(^-^) < < e^(^+^)l = 1, 

TV-foo ■' 


( 6 ) 


where V is given by the solution of the following optimization problem: 

V := inf V(x*,z) 
zeA\o 


(7) 


for V (a 


1 z) := inf It x 

T>O,<^GI3([O,T];A):0(O)=x*,0(T)=z 


^(</>) 


(see 10 


Second, we are interested in the place, where Z^ leaves the domain O. Although existing Freidlin-Wentzell 
theory cannot be applied to the epidemiological models in Examples o and |1.2| we strongly believe that the 
following holds (see 12 13 again for corresponding results for other processes): for any compact set C C 
do n (A \ O) with V < \niz(zcV{x* ,z), x G O 


lim ¥\Z^’^{t^’^) e Cl =0. 

N^oo 


In particular, if there exits a, z* G dO n (A \ O) with V{x*,z*) < V{x*,z) for all other z G dO n (A \ O), then 
for all 5 > 0, 

lim P[|Z^’"(t^’") - z*| > ,5] = 0. (8) 


1.3. Research questions 

In order to better understand the differences between deterministic and stochastic compartmental models 
in epidemiology, the theoretical results above suggest a thorough numerical analysis of the specific models of 
Examples o and |1.2| In the remainder of this article, we address the following questions. 

1. In order to analyze the deviation of the stochastic model from its Law of Large Numbers limit, we first 
have to compute the solution of the deterministic models in a reliable way in section In order to 
avoid instabilities of the solution (such as components becoming negative), we apply a non-standard 
finite difference scheme developed in [^. 

®We want to remark here that we cannot consider the exit from the domain of attraction of the disease-free equilibrium x, i.e., 

O = {z S j 4| lim Z^{t) = x}. 

t—^OO 

The reason for this is that the absorbing set A = {z £ A\zi = 0} (where 2:1 is the first coordinate of z) cannot be left both by the 
solution of the ODE and by the solution of the SDE. 

^Respectively in the time when the process is no longer attracted by the endemic equilibrium and hence the probability of 
extinction is significantly increased, cf. the LLN, equation §. 
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2. The computation of the time (and place) of exit from domain depends on the solution of the optimization 
problem Q. We use dynamic programming to solve the optimization problem numerically and hence 
deduce the results about the time of exit according to ^ in section]^ 

3. Finally, we simulate the process in section]^ It is not difficult to simulate the process in an exact way 
by using Gillespie’s stochastic simulation algorithm (SSA) (see, e.g., [l4l[T^ ). However, the simulation 
is rather slow as soon as N becomes large, and we hence adapt the t- leaping algorithm (see, e.g., [Th]). 

2. Computation of the solution of the ODE 


The large deviations theory will tell us about the deviations between the long-term behaviors of the solution 
of the SDE and of its ODE LLN limit. Hence we are interested in the long-term behavior of the deterministic 
epidemiological models described by equation (§. Consequently we need a numerical scheme to compute the 
solution of the ODE ([^ in a “reliable” way. By reliable, we mean in particular that the fixed points of the 
solution of the (discrete) numerical scheme are the same as the equilibria of the ODE (with the same local 
stability), which forces the long time bahavior of the numerical scheme to resemble that of the solution of the 
ODE. Explicit finite difference schemes for the solution of ODEs can have undesirable instabilities. In section [272| 
below, we consider the SIS model from Example |1.1| and illustrate that the solution of such a scheme can be 
oscillating around zero (and in particular become negative) if the parameters (particularly the step-size h) are 
not chosen with great care. For this reason, we apply a non-standard finite difference (NSFD) method for models 
in epidemiology which respects the qualitative behavior of the solution of the ODE (see [^). This NSFD method 
follows two basic rules for NSFD methods (for a detailed description of the NSFD method see, e.g., [I7|[I^ ). 
First, quadratic terms of the ODE are approximated in a non-local way (by mixed explicit/implicit terms). 
Second, non-trivial denominator functions are used for the discrete derivatives. 

The remainder of this section is structured as follows. In section [TH we describe the specific NSFD method 
for epidemiological models from and outline in which sense the qualitative long-term behavior of the ODE 
is respected. In section 2.2 we apply the method to the SIS model (Example [T3 and compare it to the 
corresponding explicit scheme. Finally, we apply the method to the model with vaccination (Example 1.2) 
and compute the boundary separating the domains of attraction of the stable endemic equilibrium and the 
disease-free equilibrium. 


2.1. A NSFD method for epidemiological models 

We consider deterministic models for which the ODE ([^ can be rewritten in the following form 

dZ 

— = A{Z)Z + f, Z(0) = x, (9) 

where A{Z) = (a^j(Z))i is a non-linear Metzler matrijj^and / S M'’*; usually, / represents the 
rates of birth and immigration into the different compartments (if they exist, cf. Example |1. 2 [ ). 

We follow and discretize time by fixing the time step h > 0; we write = fnh for m S N and approximate 
the solution Z of the ODE ([^ at the discrete time points tm by the solution Z of the following implicit NSFD 
scheme (i.e., « Z{tm))'- 


_ yrn ^ 

' ^(/,) ' + {i = l,...,d), Z° = x-, 

the denominator function ip is specified below and satisfies ip{h) = h o{h?). We first note that this scheme 
can equivalently be written as 


(J - iP{h)A{Z'^))Z'^+^ = -t ip{h)f, Z° = x; 


(10) 


Metzler matrix is a matrix A = (ciij)i,j=i...,d which satisfies aij > 0 for i 7 ^ j. 
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{I — 'ip{h)A{Z'^)) is an M-matri3|^and hence invertible (see [^). For the definition of the denominator function, 
we assume that all equilibria of the ODE ([^ are hyperbolic. With (^ : M —>■ M satisfying ip{h) = h + o{h?) and 
0 < i^{z) < 1 for z > 0 (e.g., ip{z) = 1 — exp(—z) or ip{z) = z/(l + z^)), V' is defined as 


ip{h) ■= 


‘pjQh) 

Q 


for Q > max 


|A|= 


2|i?e(A)| 


( 11 ) 


where the maximum is taken over the eigenvalues A of the Jacobian of the right hand side of (|^ at the equilibrium 
points. In particular, we have ipW = h + o{h?') as desired. 

The NSFD scheme (101 is elementary stable in the sense that it has exactly the same fixed points as the 
continuous ODE ([^ it approximates, with the same local stability for all values of h (see 0) In addition, it 
can be ensured that the components of the solution of the scheme remain positive by choosing h and hence 
(j){h) small enough (see again). Eor the examples we consider, the solution of the scheme always stays in the 
“natural domain” of the solution of the ODE ([^. 

2.2. The SIS model 

For the SIS model of Example |1.1[ the ODE Q becomes 

nt 


Z{t)=x+ j [j3Z{s){l — Z{s)) — ^Z{s))ds = X + f [ —l3Z{s)'^ + {/3 — j)Z{s))ds-, 

Jo Jo 


( 12 ) 


recall that in this case Z{t) denotes the number of infectious individuals at time t. This is a Riccati differential 
equation with constant coefficients which admits an explicit solution: 


f (/3-7)a;exp((;3-7)t) 

Z{t) = < (d-A+dxie^Piid-'rA)-'^) 

\ l — 0xt 


if /3 ^ 7 
else. 


(13) 


We nevertheless compute the solution of the ODE ( [T^ numerically by a standard (explicit) difference scheme 
and by the NSFD scheme introduced in section |2.1| in order to illustrate possible instabilities of the explicit 
scheme. 

Let us first consider the following explicit scheme: 


^m+1 _ 


= /3Z™(I - Z^) - yZ™, Z^ = x 


(14) 


or equivalently 

^"*+1 = Z™(I- 7 /i + /3/i)-h/3(Z™)2, Z"^=x. 

We note that if, e.g., fJ = 1/h and 7 = 2//i, Z'^ < 0 for all to > 1 which is obviously undesirable. If, e.g., 
(3 = Ijh and 7 = 4//i, the solution is oscillating around zero with Z'^ < 0 for to odd and Z'^ > 0 for to even. 
We illustrate the solution for the latter case in the left picture of Figure 2. 

Let us now apply the NSFD scheme to the model. We represent the ODE in the form (|^ by setting 
A{Z) = —f3Z+l3 — "f and define the denominator function i/i according to equation (II) (for ip{z) = I —exp(—z)). 


The solution does not suffer from the same instabilities as the solution of the explicit scheme as illustrated in the 
middle picture of Figure 2. We compare the solutions of the two schemes with the exact solution of the ODE as 


given in equation (13) (right picture of Figure 2). We want to remark here that the NSFD scheme is of course 
by no means the only possible scheme which avoids instabilities in this simple case. However, provides us 
with a theoretical basis which ensures that instabilities do not occur even for much more complicated models. 


®An M-matrix is a matrix A with aij < 0 for i ^ j which is invertible such that all entries of its inverse are non-negative. 
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Figure 2. The left picture denotes the solution of the explicit finite difference scheme (14). The 
middle picture denotes the solution of the NSFD scheme. The right picture denotes the exact 
solution of the ODE according to equation (13). Z{0) = 0.3, T = A, h = 0.1, /3 = A/h = 40, 
7 = 2//i = 20. We choose here a large step size h on purpose, to show how the two schemes 
behave for large step size. Of course, both improve if we reduce h. 


2.3. The model with vaccination 

We now consider the model with vaccination from Example |1.2[ We denote the proportions of susceptible, 
vaccinated and infectious individuals at time t by Zi{t), Z 2 {t) and Z'i{t) = 1 — Zi{t) — Z 2 {t), respectively. Note 
that we cannot reduce the dimension from three to two if we want to represent the ODE in the form (|^ (with 
appropriate Metzler matrix A). We set 


/ -/3Z3 - ^J.-r| e 7 \ 

A{Z) = I rj —aPZs 0 I and / = (/r,0,0)^ 

\ PZ 3 aPZs ) 

and choose the denominator function ij; according to equation © for (p{z) = 1 — exp(—z). We readily observe 
that the relevant assumptions for the NSFD scheme are satisfied and apply the scheme for time step h = 0.1 
to different initial values Z{Q) and a set of parameters (cf. Figure 3) such that there exists a stable disease- 
free equilibrium {x = (0.14,0.86,0)^), a stable endemic equilibrium {x* = (0.24,0.45,0.31) ^) a nd an unstable 
endemic equilibrium (i = (0.23,0.59,0.18)^; cf. the corresponding discussion in Example |l.2| ). We illustrate 
the evolution of the proportions of vaccinated and infectious individuals for different initial values in Figure 3. 
While the disease eventually dies out in the left picture, it becomes endemic in the right picture. 

The NSFD scheme also allows us to approximately compute the characteristic boundary, i.e., the boundary 
separating the domains of attractions of the two stable equilibria x and x* (in other words: the set of points 
for which the solution of the ODE converges to the unstable equilibrium x). We illustrate the characteristic 
boundary in Figure 4. 


3. Computation of V 

In order to answer the question of when the disease dies out, we have to compute the quantity V (cf. equa¬ 
tions ([^ and (0). In other words, we have to solve a control problem. We first fix the time horizon T > 0. 
We recall that we only have to consider absolutely continuous trajectories (as else It,x(4’) = 00 ) and hence we 
control the trajectory (j) via its “derivative”; for t G [0,T), s G [t,T] and an integrable function a, we write 



(j)°‘{s) := (l){s) = (j){t) + 


a(r)dr; 
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Figure 3. Evolution of the proportions of vaccinated individuals (dashed lines) and infectious 
individuals (solid lines) for initial values Z2{0)=0.5, ^3(0) = 0.05 (left picture) respectively 
2^2(0)=0.5, .^3(0) = 0.2 (right picture). T = 600, /? = 3.6, 7 = 1, 77 = 0.3, 9 = 0.02, /r = 0.03, 
a = 0 . 1 . 



Figure 4. The characteristic boundary separating the domains of attractions of a; = (0,0.86)^ 
(on the left) and x* = (0.31,0.45)^ (on the right). The unstable equilibrium x = (0.18,0.59)^ 
is on the characteristic boundary. Note the slight abuse of notation due to the reduction from 
dimension three to two. The parameters are as in Figure 3. 


for X G O, we define the set of admissible controls starting from {t,x) by 

. 4 ^ (t,a;) := |a : [t,T] —)• integrable ^ + ^ a{r)dr £ O^s G [t,T] and x + a{r)dr £ A\o'^ . 
The cost of an admissible control a starting from (t, x) is given by 

J^(t,a;, a)=^ L{(l)°‘{s),a{s))ds, 


(15) 
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and thus the value function of the optimization problem becomes (t € [0,T)) 


v^(t,x)= inf J(t,x,a) 


(16) 


and 


\T,x) = 


0 if a; ^ O 
oo else. 

The value function thus has a singularity at the terminal time T which is due to the constraint (^(T) ^ O. For 
the models we have in mind (cf. Examples o and |1.2[), i t is possible to move on dO n (A \ O) by following the 
solution of the ODE ([^ (which is free of cost, cf. ([^)p°| Eor this reason, we have u^(t, x) < x) for T < T 


and therefore 


lim 11^(0, X*) = V. 
T —¥OC) 


We want to remark here that any admissible control a with </>“ (T) ^ (A \ O) n dO (cf. (§ again) cannot be 
optimal. Therefore, it would be equivalent to require the weaker restriction x + a(r)dr € A for s € [0,T] in 
equation (15). 

In the following, we describe a discrete dynamic programming algorithm for the numerical approximation 
of V (section |3.l[ ) and apply this algorithm to the models from Example (section |3.2[ ) and Example |l.2| 
(section |3.3[ ). 

3.1. Approximation by discrete-time dynamic programming 

In order to numerically compute the solution of the optimization problem ( [l^ , we apply discrete-time 
dynamic programming. To this end, we discretize time and consider n -|- 1 (n € IN) equidistant time points 

T 

0 = to < ■ ■ ■ < tn = T, i.e., =tm = mAt for At = —. 


For j = 0, ..., n — 1 let q; = (5(tm))m=j,....n-i, where a{tm) & We consider piecewise constant trajectories 
and dehne recursively for x S O, 

^°‘{tj) = ^{tj) = X, ^{tm+l) = Htm) + a{'tm)At {m = j, . . . , U - 1). 

We define the set of admissible controls by 

A^'^{tj,x) ;= |a (f>°‘{tjn) & dVm = j,..., n — 1 and ^“(t„) G (d \ O) n do|. 

The cost of an admissible control a G A^’"‘{tj,x) is given by 

n—1 

J^’”(t^-,X, a) At ^ L{(t)°‘{tm),a{tm)) 


and the value function is given by 

v'^’'^{tj,x) := inf j'^'^{tj,x,a). (17) 

,x) 


^^For the SIS model, we have dO fl (A \ O) = {0}. For the model with vaccination, we have 

dO n (A \ O) = G A| lim Z^{t) = x} 

t—^oo 

(recall that x is the unstable endemic equilibrium). 
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The constraint implies that 


= At inf L{x,a). (18) 

Q!GR‘^:a:+Q:AtG{A\0)n50 

Furthermore, for j < n — 1, the dynamic programming principle implies the following Bellman equation: 
v^'"'{tj,x) = inf |l(x, a)At + a; + Q;At)|. 

aGlEl'^:x+QAtGA J 

Hence the solution of the optimization problem can be computed backward in time. 

We need next to discretize space. Note also that in the SIS model, the functional L is explicitly known, while 
this is not the case in the vaccination model. For that reason we now distinguish between the two cases. 

3.2. The SIS model 

We discretize space via a grid of the set A = [0,1] with fineness Ax = lln^ for some new integer h. Hence 
the grid points are given by 

x'' := iAx for i = 0,..., h. 

Given a function / which is only defined on the grid points, we define values of the function for arbitrary 
X £ Ahy its linear interpolation /int. In other words, we obtain for a: S H, 


n 

fint{x) := l{,rg[,^._,,.+i)}|/(a:*)(a;*+^ - a:) + /(a:*+^)(a: - a:*)|. 


i=0 


We obtain the following (backward in time) algorithm for the solution of the optimization problem (17). 

1. For X on the grid. 


a;) = At inf L{x,a). 

a^^^:x-\-aAt^{A\0)r]dO 

We denote the solution of this minimization problem by a*{t„-i) = a*{tn-i,x). 

2. For m < n — 1, X on the grid, 

i;^’”’"(f„,x) = inf (L{x,a)At + v'[^^’"^{tm+i,x + aAt)\. 

We denote the solution of this minimization problem by a*{tm) = (x*{tm,x). 

Note that this algorithm involves the solution of n x h minimization problems. We recover an approximation 
(j)* = {(j)*{tm))m=Q,...,n of the Optimal trajectory starting from the stable equilibrium x* , by the following 
procedure (including another linear interpolation). 

3. For m = 0,..., n — 1, 

0*(to) = a:*, = 4>*itm) + (j)* {tm))At. 


It is well-known that under appropriate assumptions, the solution of this procedure converges to the solution 


of the optimization problem (16) (see, e.g., [^[^); the optimal coupling of time and space discretization is 
At = Ax. We want to remark here that standard assumptions do not apply to our model. For instance, the 
function L is not Lipschitz continuous, even L(x,y) —> oo for x —)■ dA is possibleFurthermore, the terminal 


^^This is e.g. the case for y <0 and a; —^ 0. Note however that the action required to approach the boundary, e.g., at constant 
speed one is finite and we have L(t — s, —l)ds —>■ 0 as t —> 0. 
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Figure 5. Optimal trajectory (j)* for T = 5 (left picture), T = 20 (middle picture) and T = 60 
(right picture); At = Aa: = /3 = 1.5 and 7=1. 

constraint results in a singularity of the value function at terminal time T. We are not aware of any analytical 
results including our set-up. We nevertheless have reasons to believe that 

V « a;*) for large enough T. 

It is not in the scope of this article to provide a thorough analytical proof of the convergence of the algorithm. 
We readily observe that for the SIS model, the Legendre-Fenchel transform is given explicitly by 

L{x,y) = j/log0- /3x(l - x){e - 1) - ix{\ - l) , 


where 


e := 


y + + 4/l7a;2(l — x) 

2j3x{l — x) 


Furthermore, the boundary consists of the single point {0}; hence equation (18) becomes 

v'^’'^{tn-l,x) = AtL{x,—x/At). 

We consider the following parameters: /3 = 1.5, 7 = 1; in particular, we have i?o = 1.5 > 1 and the endemic 
equilibrium x* = 1/3 is stable. We apply the dynamic programming algorithm described in section 3.1 with 
At = Ax = 1/100 and T = 5,10,20,40,60 (i.e., n = 100 and n = 500,1000,2000,4000,6000). We obtain 


V: 


®’"’"(to,a:*) = 0.0953, a:*) = 0.0757, uf„°’”’"(to, a:*) = 0.0705, 


int 

X*) = 0.0702, C’"’”(to, X*) = 0.0702. 


int 
, 60 ,r 
int 


We hence deduce that V « 0.0702 and thus « exp(0.0702iV) for large N. We illustrate the solution of 

the optimization problem in Figure 5 for T = 5, 20, 60. We observe that it is the cheapest to move slowly near 
the stable equilibria x* and the unstable equilibrium x = 0. 


3.3. The model with vaccination 


Unlike in section 3.2 the place of exit from O is a priori not clear for the model with vaccination (cf. Exam¬ 
ple 1.2). It seems plausible that it is cheapest to exit O at the unstable endemic equilibrium x. Indeed, if T 
is chosen large enough, an optimal (or near optimal) trajectory (j) can move on the boundary towards x after 
the exit time without generating any extra cost (cf. (§); in this case we have limt_>oo />(t) = x. Hence, we can 
assume that x is indeed the place of exit for the numerical computation of V. Therefore, equation (18) becomes 


v'^’"'{tn-l,x) = AtL{x, (x — x)/At). 
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However, the Legendre-Fenchel transform L cannot be computed explicitly. Here we use the following defi¬ 
nition of L, which is equivalent to Q (see [To]): 


where for /i S 


L{x,y) := 


inf^GBx.B if7^ 

00 if not, 


{Pj (x) - H + H log (^)}; 


and 


B. 


.^y := S > 0 only if I3j{x) > 0 and y = ^ 


Since we seek to minimize 

V := inf [ cj)'(t))dt, 

s.t. (/>(0)—x* ,(/>(T)—:r Jq 

our discretized dynamic programming equation reads: for each x € Q (the space grid), 

v{tk,x) = inf \l[x,— 
x'Gy ^ \ L 

= inf^ ^£{fL{x,x'),x)At + v{tk+i,x')'j, 
where jl{x,x') is the argument which realizes the following minimum 


— X 

At 


At- 


;(4+i,x')| 


inf i{fj,,x). 

ti, x+A J2j IJ-jhj=x' 

Note that here we avoid to consider interpolation between points on the grid. We have implemented this 
algorithm on a grid whose one of the axis is the line passing through x* and x, with the following values of the 
parameters : /3 = 3.6, 9 — 0.02, y = 0.03, 17 = 0.3, a = 0.1 and 7 = 1. 

The approximation of V which we have found is 0.3891, with At = 0.05, and the Axi = Ax 2 = 0.005. 
This means that ] « exp(0.3891A) for large N. Note that if At is not small enough compared to Ax, 

since we restrict the above minimization to cc' € G, we cannot allow small speed, and we get larger values of 
V. We expect that more accurate results could be obtained, by doing the minimization over x' on a Hner grid, 
which necessitates to interpolate v(tk+i,-) between points on the grid. We have not been able to implement 
this refinement, by lack of time. This will be done in the future. It should improve the numerical estimate of 
V. The above quantity should be considered as an upper bound of the true value. 

4. Simulation of 

In this section, we describe known exact and approximate simulation methods for the simulation of the 
solution of equation Q. To this end, we fix A; with the notations ^ and aj{z) = Nj3j{z), equation Q 
can be rewritten as 

Z^{t) = x + '^iZjPj f aj{Z^is))ds\. 

3=1 ^ 


(19) 
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Note that our ultimate goal would be to simulate the value of which, from the above Large Deviations 
results, is huge, as soon as N is large. This may be seen as impossible. However, this is similar to the problem 
of rare event simulation, for which there exist specific techniques, which we intend to implement in the future. 
See also the discussion below in section |4.7.2[ In the framework of this paper, we present some preliminary 
numerical work, which aims only at simulating the stochastic process over a given finite time horizon. 


In th e fol lowing, we first introduce the so-called stochastic simulation algorithm (SSA) proposed in 14 15 
(section 4.1). The SSA allows for the exact simulation of the process . Unfortunately, this algorithm is 


rather slow in some circumstances, typically when there are too many jumps of the Poisson processes on a 
fixed time interval. Therefore, we discuss faster approximate simulation methods in the subsequent sections: in 

we discuss 


4.2 we introduce the explicit r-leaping algorithm (see 16 ) in the form of in section 4.3 


section 

modifications of the r-leaping algorithm: implicit and mid-point r-leaping (section 4.3.1) and the algorithm 
from 22 which ensures that the simulated process remains in its domain, e.g., by preventing the components 
of the process from becoming negative (section [4.3.2[ ). We apply the SSA and the modified r-leaping algorithm 
from section |4.3.1| to the specific models from Examples o and |1.2| in sections [4(4| and [451 respectiv ely. We 
conclude by addressing open research questions concerning the simulation of the process Z^ in section 4.7 


4.1. The stochastic simulation algorithm 


A mathematically exact procedure for simulating the evolution of the process Z^ given by equation (19) is 


the stochastic simulation algorithm (SSA) first proposed in 14,15 . The SSA changes the value of the solution 


at each jump time of one of the Poisson processes. The simplest implementation of the SSA is the so-called 
direct method: 

Whenever t <T: 

1. While in state z at time t, evaluate all functions aj{z) and their sum ao(z) := 

2. Generate the time increment r as the realization of an exponential random variable with parameter 
ao{z). 

3. Choose an index j G {1,..., fc} according to the probability distribution { j = 1,..., fc). 

4. Update 

i. t t + T, 

ii. z <— z + Vj. 

5. Record (t, z). Return to 1. if t -I- r < T; else stop. 

Carrying out step 3. is mathematically straightforward: we draw one realization r of a uniform random 
variable on [0,1]; then, j is the smallest positive integer for which 


'^Uiiz) > rao(z). 

i=l 

The SSA is mathematically exact. However, the task of explicitly simulating each jump makes the SSA too 
slow for practical implementation in some circumstances. The reason is that for large N, there are too many 
jumps on any fixed time interval. 

4.2. The explicit Poisson r-leaping algorithm 

A faster but approximate stochastic simulation procedure is the explicit Poisson r-leaping algorithm (see [16]). 
The basic idea of this procedure is to advance the system by a preselected time increment r (in contrast to the 
randomly generated time increment r in the SSA). On the one hand, r should be chosen large enough in order 
to ensure that “many” jumps occur during this length of time (as else, the savings in computation time are not 
significant); on the other hand, it should be small enough such that the values of the functions aj are unlikely 
to change “significantly” (as else, the procedure is not accurate enough). The latter restriction is called the leap 


condition. A strategy due to 21 for satisfying this condition is to require that the changes of the functions Oj 
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during a leap r are “likely” to be bounded by eao(z); here, e (0 < e 1) is the error control parameter. In 
other words, 


\aj{Z" {t + t)) — aj{Z" (t))| < eao{Z" (t)) for j = 1..., A: with high probability. 


We estimate the largest value of r that meets this particular requirement as follows. First, we compute 

d 


fjj' (^) := for ■ • •: 


2 = 1 


dzi 


( 20 ) 


here Vij denotes the i-th component of the vector Vj. Subsequently, we calculate 

k k 

:= and a]{z) := ^ f]y{x)aj{z) for j = 1,..., fc 

i'=i i'=i 


and set 


= r{z) := 


eao(z) 




i=.i....fc \^ij{z)\ aUz) 


( 21 ) 


( 22 ) 


Note that for a standard Poisson process P, 


aj{Z^is))ds^ - a,{Z^{s))d. 


s]=P 


aj{Z^ {s))ds^ 
Piaj{Z^(t))T) = P{a^{z)T) =: Pj 


and thus 

k k 

Z^it + r) « Z^{t) p'^PjVj = z P'^PjVj. 
i=i 

This procedure for the selection of r yields the following r-leaping algorithm. 

Whenever t <T: 

1. While in state 2 at time t, evaluate the functions Oj and their sum ao(^) = X)hi 

2. Compute r via the formulas (20l-(22). 

3. Fix a small integer n (typically, e.g., n = 10) and a moderate integer n (typically, e.g, n = 100). 

i. If T < reject r and execute n iterations of the SSA. Then, return to step 1. 

ii. If T > proceed to step 4. 

4. For j = 1, k, generate pj as a Poisson random variable with mean aj{z)T. 

5. Update 

i. t ^ t + T, 

ii. z ^ z + Y!"^^^pjVj. 

6. Record {t, z). Return to 1. if t + r < T; else stop. 

Let us shortly comment on this algorithm, pj represents an approximation of the number of jumps of type j 
during the time interval [t, t+ r[. We remark again that the approximation by Pj in step 4. is justified as long as 
aj{z) remains more or less constant over the subsequent time interval of length r. The reason for step 3. i. is the 
following, is the mean time step to the next jump when applying the SSA. If r is not significantly larger 

than the computation of r (via equations (20)-(22)) is no more efficient then applying the SSA directly 

(which is also of course more accurate); as this property is likely to persist for a while, it is more appropriate 
to apply the SSA for some time before returning to the r-selection procedure again. 
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The explicit r-leaping algorithm above has been shown to provide accurate sim ulation re sult s; furthermore, 
it is substantially faster than the SSA for many “not-too-stiff” system^ (see 


22 


see also 


24 


for theoretical 


evidence for this). We illustrate the gain of computation time by numerical experiments in sections 4.4 and 4.5 
below for a modification of this algorithm (cf. section 4.3.2). 


4.3. Modifications of the r-leaping algorithm 

There exist several variants of the above described explicit r-leaping algorithm. We now describe some of 
these. 


4.3.1. Implicit and mid-point r-leaping algorithms 

One may want to freeze the rates not at the value of the solution at the beginning of the discretization 
interval, but rather at an approximate value of the position at the end of that interval, or else at the mid-point. 
The “implicit” r-leaping algorithm is obtained by replacing aj{Z^ {t))T with 

-br6(Z^(t)))r, 


and the mid-point r-leaping algorithm consists in replacing the same quantities by 

aj{Z^{t)+Tb{Z^{t))/2)T. 


In 24 it is proved that the error in the mid-point r-leaping algorithm is less than the one produced by the 


explicit r-leaping algorithm. 

4.3.2. A modified r-leaping algorithm 

As most numerical schemes, the r-leaping algorithm has the disadvantage that the approximate solution 


might very well jump over the boundaries which the solution of the stochastic equation (19) can never cross. 
More precisely, the solution of 0 satisfies certain constraints. For example, in the case of constant total 
population, we have that Z^ {t) > 0, 1 < i < d, and Z^ {t) := 1 — ^ 0- We would like that the 

approximate solution satisfies these constraints as well. However, this might not be the case with a standard 
r-leaping algorithm. In particular, if for some 0 < * < d, Z[^{t) is very small, then some of the Poisson driving 
processes may push it over to a negative value during the next time interval of length r. Because of that, at 
the beginning t of each time step, we define 


L,:= 


0 <i<d, Uijaj{Z^ {t))<.0 


Ziit), 


d 

where voj := — Vij, 

i=l 


and choose a critical value Uc- The critical directions will be chosen by comparing Lj to Uc- 

The modified Poisson r-leaping procedure requires that a critical Poisson process cannot jump more than once 
in a time interval of length r. That makes it impossible for a critical Poisson process to drive any component of 
Z^ to a negative value. The theoretical justification for each step in this modified r-leaping procedure is given 


in 


22 


Whenever t <T: 


1. While in state z at time t, evaluate the functions aj{z) and their sum ao(z) = 'Z!j=i 

2. Identify the currently critical Poisson processes as those Pj for which 


aj(z) > 0 and Lj < Uc- 

^^I.e., systems for which the difference between the smallest and the largest jump rates is not too large (see, e.g., 221). 
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3. Modify the r-selection procedure from equations (20)-(22| by only taking into account those jump 


4. 


5. 


6 . 


directions j' corresponding to the non-critical Poisson processes. Denote the resulting time step by r'. 
If there is no non-critical Poisson process, set r' = oo. 

Fix a small integer n (typically, e.g.,n = 10) and moderate integer n (typically, e.g, h = 100). 

i. If r' < reject r' and execute h iterations of the SSA. Then, return to step 1. 

ii. If t' > proceed to step 5. 

Compute the sum a(j(z) of the functions aj(z) of the critical Poisson processes. Generate r" as a 
realization of the exponential random variable with parameter aQ{z). 

i. If t' < r", set t = t' . For all the critical Poisson processes set pj = 0. For all the non-critical 
Poisson processes generate pj as a realization of the Poisson random variable with mean aj{z)T. 

ii. If t" < t', set t = t" . Generate an integer according to the probability distribution aj( 2 ;)/ag(z) 
(i.e., j runs over the index values of the critical Poisson processes). Set pj^ = 1, and for all the 
other critical Poisson process, set pj = 0. For all the non-critical Poisson processes Pj, generate pj 


as a realization of the Poisson random variable with mean aj{z)i 


7. Update 


i. t t + T, 

ii. z ^ z + J^UiPj'^P 

If z is not in its domain, replace t' ^ t'/ 2, and return to step 6. 
Record (t,z). Return to 1 if t -I- t < P; else stop. 


Remark 4.1. It may seem strange to use a uniform critical value ric, independent of the value of the rate 
of each Poisson process. It would be natural to compare each Lj with the corresponding aj{Z^ {t))T, which 
however would require to compute a first value of r. We intend to search for improvements in that direction in 
the future. 


Remark 4.2. Although the modified r-leaping algorithm is designed so as to reduce the chance of a negative 
value of any of the it), it cannot completely rule them out. This is the reason for step 8. of the algorithm. 
We suspect that this step might introduce a bias in the simulation, in that it tends to delete simulations 
with exceptionally large values of some of the Poisson process increments, and also to reduce the chance of 
approaching or hitting any of the boundaries, which is specially annoying for our purposes. 

We intend to check two other possible approaches. The first one consists in projecting the computed point, if 
it is outside the domain. The second consists in replacing r by r/2, but without resimulating a new independent 
value of the Poisson process increment, but rather by simulating the new increments according to the conditional 
law of Pj{t + T') — Pj{t), given the previously simulated value of Pj{t + T) — Pj{t), which, if Pj{t + T) — Pj{t) = k, 
follows a binomial law with parameters k and r'/r. 


4.4. The SIS model 


In this section, we apply the simulation algorithms from sections 4.1 and 4.3.2 to the SIS model (Example |l.l[ ). 
We consider an endemic situation {(3 = 1.5, 7 = 1, hence Rq = 1.5 > 1 and x* = 1/3) and compare the 
performance of the SSA to the modified r-leaping algorithm. If the population size is relatively small, N = 2000, 
we observe that the computation times are similar for both algorithms (even smaller for the SSA). A closer look 
at the realized trajectory of the r-leaping procedure confirms that almost all jumps are of size 1/2000. This 
suggests that in most cases, the selected time step r' is rejected and the SSA is applied (cf. steps 3. and 4. i. of 
the algorithm). Hence, there is no significant gain in computation time. If we let N = 20000, the number of 
jumps (and hence the computation time) for the SSA increases approximately by a factor ten. The r-leaping 
algorithm is significantly faster which suggests that the time step r' is now rarely refused. While the SSA 
becomes inefficient for very large time horizons (e.g., T = 600) and population sizes (e.g., N = 200000), the 
process can still be simulated efficiently by the r-leaping algorithm. 
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Figure 7. Realized trajectories of the modified r-leaping procedure for N = 2000 (left picture) 
and N = 20000 (right picture). The parameters are as in Figure 6. 


We illustrate realized trajectories of the SSA and the modified r-leaping algorithm in Figures 6 and 7, 
respectively (left pictures: N = 2000, right pictures: N = 20000). We observe that they approach the solution 
of the ODE as N increases for both algorithms as predicted by the LLN (cf. (|^ ; see i) 


4.5. The model with vaccination 


We now apply the same algorithms to the model with vaccination (Example |1.2[ ). We consider the bistable 
situation from section 2.3 (cf. Figure 3 for the parameters) and observe the same qualitative properties concerning 
computation times as for the SIS model in section [T4| In Figures 8 and 9, we illustrate the realized trajectories 
for the SSA and the modified r-leaping algorithm, respectively (left pictures: N = 2000, right pictures: N = 
20000). For simplicity, we only display the proportions of infectious individuals. As in section 4.4 we observe 
that the trajectories approach the solution of the ODE as N 


oo. 


It is tempting to compare those figures, in order to obtain the error induced by the tau-leaping scheme, 
compared to the exact SSA scheme. However, this would require that the two schemes are computed using the 
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Figure 8 . Realized trajectories for the proportion of infectious individuals simulated by the 
SSA (left picture: N = 2000, right picture: N = 20000). ^2(0) = 0.2, ^3(0) = 0.7, T = 50. All 
other parameters are as in Figure 3. 




Figure 9. Realized trajectories for the proportion of infectious individuals simulated by the 
T- leaping algorithm (left picture: N = 2000, right picture: N = 20000). All parameters are as 
in Figure 8. 


same realization of the Poisson processes. This is not what we have done here. We intend to do this in a near 
future, in order to confront simulations with the predictions of the theory concerning the error induced by the 
tau-leaping algorithm, see e.g. 


23 , 24 


4.6. Computation time 

In the next table we present the computation times (in seconds) of the SSA and modified tau-leaping for 
different population sizes in the SIS model. All programs were written and executed in Matlab. We remark 
that from N = 20000 on, the computation time of the tau-leaping scheme is spectacularly smaller than that of 
the SSA scheme. We also note that the computation time of the tau-leaping scheme seems to be a decreasing 
function of N. The reason is probably that the larger N, the less frequently the scheme uses SSA steps. 
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SSA 

Population size N 

Time horizon T 

Computation time 


2000 

50 

2.8 


20000 

50 

24.66 


200000 

50 

225.44 

Modified tau-leaping 





2000 

50 

2.3 


20000 

50 

3.02 


200000 

50 

1.05 


The next table presents the computation times of the SSA and modified tan-leaping for different population 
sizes in the SIV model. The same comments as for the SIS model apply. 


SSA 

Population size N 

Time horizon T 

Computation time 


2000 

50 

4.51 


20000 

50 

36.4 


200000 

50 

352.64 

Modified tau-leaping 





2000 

50 

2.95 


20000 

50 

1.04 


200000 

50 

0.67 


4.7. Open questions 

We conclude by describing two further open questions concerning the simulation of the process . We 
discuss briefly the simulation of models where some jump types are significantly more frequent than others 
(section 4.7.1). Finally, we discuss the simulation of rare events with the help of large deviations theory 


(section 4.7.2). 


4.7.1. Models with different time scales 

The simulation algorithms described in sections |4.1||4.3| tend to become inefficient or inaccurate if some of 


the jump rates are much larger than others (cf. also Footnote 121. In epidemiology, this is typically the case for 
so-called host-vector diseases (e.g., malaria), where the disease is transmitted from vectors (mosquitoes) to hosts 
(humans) and vice versa. In the particular case of malaria, the jumps concerning only the mosquito population 
have a much shorter time scale than those also concerning the human population (see, e.g., the model in [^). 
It is hence time consuming to simulate a significant number of jumps within the human population. 

In order to overcome this problem, we simulate only the jumps with small rates via r-leaping. Within the 
interval [t,t + r[^ the jumps with large rates are approximated by the solution of the ODE (cf. section 
We intend to analyze the efficiency and accuracy of this procedure. In particular, we are interested in whether 


the large deviations behavior (cf. section 1.2) is changed by this partial ODE approximation signihcantly. 


4.7.2. Simulation of rare events 

We are particularly interested in sampling the probability of the event that the process exits the domain 


of attraction of a stable equilibrium before a given time T (cf. section 1.2). If iV is large and T is relatively 
small, the probability of this event is typically extremely small. Applying a naive Monte Carlo (MC) approach 
is therefore inefficient and we rely on rare events simulation techniques. In 25 , the application of Importance 


Sampling is explained in the context of diffusion processes. We intend to adapt this method in the following 
way. Via a change of measure, we transfer the process Z^ to a jump process Z (see, e.g., for a Girsanov 
theorem for processes of the type of Z^) for which the event of interest is much more probable. Then, we apply 
MC simulation to the process Z and associate a weight with each realization (hence the average in the MC is 


13 


Here, r is the time step chosen by the r-selection procedure applied to only the small rates. 

^'^The initial value 2 : for this computation ca n eithe r be chosen explicitly, i.e., 2 : = Z^{t) or implicitly/at the mid-point of the 
interval (cf. the somehow related issue in section 


4.3.Id 
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replaced by a weighted average). In order to choose the process Z appropriately (i.e., such that the probability 
of exit is “maximized”), we use the theory of large deviations: the “cheapest” way to exit the domain is given 
by the trajectory (cf. sections 1.2 and[^; hence, we change the rates of the process in such a way that the 
corresponding rate function /(<))*) becomes small (ideally, 0* is approximately the solution of the ODE for the 
modified rates). 
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